home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / nalgfa.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  13.4 KB  |  458 lines

  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
  9. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10.  
  11. (in-package "MAXIMA")
  12. (macsyma-module nalgfa)
  13.  
  14. (declare-top(special alg-num vlist *var *nosplitf *algvar 
  15.           *denom *num *ans algfac* $nalgfac ALPHA)
  16.      (genprefix nalg))
  17.  
  18. (load-macsyma-macros rzmac ratmac)
  19.  
  20. (setq alg-num 0)
  21.  
  22. (defun new-alg nil
  23.        (newsym (concat '$alg  (setq alg-num (f1+ alg-num)))))
  24.  
  25.  
  26. (defun psqfrp (p var)
  27.        (zerop (pdegree (pgcd p (pderivative p var)) var)))
  28.  
  29. (defun pgsubst (val var p)    ;;generalized psubst substitutes any
  30.        (cond ((pcoefp p) p)    ;;expression for any var in p
  31.          ((eq var (car p))
  32.           (cond ((pzerop val)
  33.              (pterm (cdr p) 0))
  34.             ((do ((ld (cadr p) (car a))
  35.               (a (cdddr p) (cddr a))
  36.               (ans (caddr p)
  37.                    (pplus
  38.                 (ptimes ans
  39.                     (pexpt val
  40.                            (f- ld (car a))))
  41.                 (cadr a))))
  42.              ((null a) (ptimes ans (pexpt val ld)))))))
  43.          ((pointergp var (car p)) p)
  44.          ((do ((a (cdddr p) (cddr a))
  45.            (ans (ptimes (list (car p) (cadr p) 1)
  46.                 (pgsubst val var (caddr p)))
  47.             (pplus ans
  48.                    (ptimes (list (car p) (car a) 1)
  49.                        (pgsubst val var (cadr a))))))
  50.           ((null a) ans)))))
  51.  
  52. (defun pvsubst (nvar ovar p)
  53.        (cond ((or (pcoefp p) (pointergp ovar (car p))) p)
  54.          ((eq ovar (car p))
  55.           (cons nvar (cdr p)))
  56.          (t (pgsubst (make-poly nvar) ovar p))))
  57.  
  58. (defun ordervar (var l)
  59.   (let ((mvar (lmainvar l)))
  60.     (cond ((null mvar) l)
  61.       ((null (pointergp mvar var)) (cons var l))
  62.       ((let ((newvar (gensym)))
  63.         (setq genvar (append genvar (list newvar)))
  64.         (setplist newvar (symbol-plist var))
  65.         (set newvar (f1+ (symbol-value mvar)))
  66.         (cons newvar
  67.               (mapcar #'(lambda (p) (pvsubst newvar var p))
  68.                   l)))))))
  69.  
  70. (defun lmainvar (l)        ;;main var of list of poly's
  71.        (do ((l l (cdr l))
  72.         (v))
  73.        ((null l) v)
  74.        (cond ((pcoefp (car l)))
  75.          ((null v) (setq v (caar l)))
  76.          ((pointergp (caar l) v)
  77.           (setq v (caar l))))))
  78.  
  79. (defun presult (p1 p2 var)        ;;change call in algsys?
  80.     (let ((genvar genvar))
  81.        (setq var (ordervar var (list p1 p2))
  82.          p1 (cadr var)
  83.          p2 (caddr var)
  84.          var (car var))
  85.        (cond ((zerop (pdegree p1 var))
  86.           (cond ((zerop (pdegree p2 var)) 1)
  87.             ((pexpt p1 (cadr p2)))))
  88.          ((zerop (pdegree p2 var))
  89.           (pexpt p2 (cadr p1)))
  90.          ((resultant p1 p2)))))
  91.  
  92. (defun pcoefvec (p)
  93.    (cond ((pcoefp p) (list p))
  94.      ((do ((l)
  95.            (i (cadr p) (f1- i))
  96.            (p (cdr p)))
  97.           ((signp l i) (nreverse l))
  98.         (push (cond ((and p (f= (car p) i))
  99.              (prog1 (cadr p) (setq p (cddr p))))
  100.             ( 0 ))
  101.           l)))))
  102.  
  103. (defun algtrace1 (bvec tvec)
  104.        (do ((i (f- (length tvec) (length bvec)) (f1- i)))
  105.        ((f= i 0) (algtrace* bvec tvec))
  106.        (setq bvec (cons 0 bvec))))
  107.  
  108. (defun algtrace* (bvec tvec)
  109.        (do ((b bvec (cdr b))
  110.         (tr (car (last bvec))
  111.         (pplus tr (car (last b)))))
  112.        ((null (cdr b)) tr)
  113.        (or (pzerop (car b))
  114.            (do ((l (cdr b) (cdr l))
  115.             (tv tvec (cdr tv)))
  116.            ((null l))
  117.            (rplaca l (pdifference (car l)
  118.                       (ptimes (car b) (car tv))))))))
  119.  
  120. (defun algtrace (r p)
  121.        (cond ((eq (caar r) (car p))
  122.           (ratreduce (algtrace1 (pcoefvec (car r))
  123.                     (cdr (pcoefvec p)))
  124.              (cdr r)))
  125.          ((ratreduce (pctimes (cadr p) (car r))
  126.              (cdr r)))))
  127.  
  128. (DEFMFUN $algtrace (r p var)
  129.   (let ((varlist (list var))            
  130.     (genvar nil))
  131.     (rdis* (algtrace (rform r) (car (rform p))))))
  132.  
  133. (defun good-form (l)             ;;bad -> good
  134.        (do ((l l (cddr l))
  135.         (ans))
  136.        ((null l) (nreverse ans))
  137.        (push (cons (cadr l) (car l)) ans)))
  138.  
  139. (defun bad-form (l)            ;;good -> bad
  140.        (mapcar #'(lambda (q) (list (cdr q) (car q))) l))
  141.  
  142. (DEFMFUN $algfac n
  143.        (cond ((f= n 3) (apply '$pfactoralg (listify n)))
  144.          (t (let ((varlist))
  145.           (cond ((f= n 2)
  146.              (newvar (arg 2))
  147.              (cond ((alike1 (arg 2) (car varlist))
  148.                 ($pfactoralg (arg 1) nil (arg 2)))
  149.                    (t ($pfactoralg (arg 1)
  150.                            (arg 2)
  151.                            (car (last varlist))))))
  152.             (t (newvar (arg 1))
  153.                (setq varlist
  154.                  (mapcan #'(lambda (q) (cond ((algpget q)
  155.                                   (list q))
  156.                                  (t nil)))
  157.                      varlist))
  158.                (cond ((f= (length varlist) 1)
  159.                   ($pfactoralg (arg 1)
  160.                            nil
  161.                            (car varlist)))
  162.                  ((f> (length varlist) 1)
  163.                   (merror "too many algebraics"))
  164.                  (t (merror "no algebraics")))))))))
  165.  
  166. (DEFMFUN $pfactoralg (f p alg)
  167.        (let ((varlist (list alg))
  168.          (genvar) (vlist) (tellratlist) ($ratfac)            
  169.          ($gcd '$alGebraic) 
  170.          ($algebraic) ($ratalgdenom t)
  171.          (*denom 1) (*num 1) (*ans))
  172.         (cond ((and (null p) (radfunp alg t)) (newvar (cadr alg)))
  173.           (t (newvar p)))
  174.         (newvar1 f)
  175.         (cond ((null vlist) (merror "attempt to factor a constant")))
  176.         (setq varlist (nconc varlist (sortgreat vlist)))
  177.         (cond (p (setq p (cadr (ratrep* p)))
  178.              (push (cons alg (mapcar (function pdis) (cdr p)))
  179.                tellratlist))
  180.           (t (setq p (algpget alg))
  181.              (setq p (pdifference
  182.                   (pexpt (cadr (ratrep* alg)) (car p))
  183.                   (cadr p)))))
  184.         (setq $algebraic t)
  185.         (setq f (cadr (ratrep* f)))
  186.         (setq f (pfactoralg1 f p 0))
  187.         (cons '(mtimes)
  188.           (cons (rdis (ratreduce *num *denom))
  189.             (mapcar 'pdis f)))))
  190.  
  191. (declare-top(special adn*)) ;also redefine fact5 to call nalgfac correctly
  192.  
  193. (defun nalgfac (p mp)
  194.        (let ((*num 1) (*denom 1) (*ans) (algfac*) ($nalgfac)
  195.          ($algebraic t) ($gcd '$algebraic))
  196.         (setq p (pfactoralg1 p mp 0))
  197.         (setq adn* (times adn* *denom))
  198.         (cond ((equal *num 1) p)
  199.           (t (cons *num p)))))
  200.  
  201. (setq *nosplitf t)
  202.  
  203. (defun pfactoralg1 (f p ind)
  204.   (cond ((pcoefp f) (setq *num (ptimeschk f *num)) *ans)
  205.     ((f= (cadr f) 1) (setq f (pshift f (car p) ind))
  206.             (push (algnormal f) *ans)
  207.             (setq f (rquotient f (car *ans))
  208.                   *denom (ptimeschk (cdr f) *denom)
  209.                   *num (ptimeschk (car f) *num))
  210.             *ans)
  211.     ((equal (cdr f) (cdr p))
  212.      (push (pdifference (make-poly (car f)) (make-poly (car p))) *ans)
  213.      (setq f (rquotient f (car *ans))
  214.            *denom (ptimeschk (cdr f) *denom))
  215.      (pfactoralg1 (car f) p ind))
  216.     ((zerop (pdegree f (car p)))
  217.      (mapc #'(lambda (q) (pfactoralg1 (pshift q (car p) -1) p (f1+ ind)))
  218.            (let (($algebraic nil)
  219.              ($gcd (car *gcdl*)))
  220.             (pfactor1 f)))
  221.      *ans)
  222.     (t (do ((l (let (($algebraic nil)
  223.              ($gcd (car *gcdl*)))
  224.             (pfactor (algnorm f p)))
  225.            (cddr l))
  226.         (polys)
  227.         (temp)
  228.         (alg (car p)))
  229.            ((null l)
  230.         (setq *num (ptimeschk f *num))
  231.         (mapc #'(lambda (q) (pfactoralg1
  232.                      (pshift q alg -1) p (f1+ ind)))
  233.               polys)
  234.         *ans)
  235.            (cond ((pcoefp (car l)) nil)
  236.              (t (setq temp (cond ((null (cddr l)) f)
  237.                      (t (pgcd f (car l)))))
  238.             (cond ((f= (cadr temp) 1)
  239.                    (setq temp (algnormal temp))
  240.                    (push (pshift temp alg ind) *ans))
  241.                   ((f= (cadr l) 1)
  242.                    (setq temp (algnormal temp))
  243.                    (push (pshift temp alg ind) *ans)
  244.                    (or *nosplitf
  245.                    (setq *nosplitf
  246.                      (list (car l) temp ind))))
  247.                   (t (push temp polys)))
  248.             (setq f (rquotient f temp)
  249.                   *denom (ptimeschk (cdr f) *denom)
  250.                   f (car f)))) ))))
  251.  
  252. (defun pshift (f alg c)
  253.        (if (= c 0) f
  254.        (pgsubst (pplus (make-poly (car f)) (pctimes c (make-poly alg)))
  255.             (car f) f)))
  256.  
  257.  
  258.  
  259. (DEFMFUN $splitfield (p var)
  260.        (let ((varlist)
  261.          (genvar)
  262.          (genpairs)
  263.          (*algvar)
  264.          ($gcd '$ALGEBRAIC))
  265.      (newsym var)
  266.      (setq *algvar (caar (new-alg)))
  267.      (setq p (psplit-field (cadr (ratf p))))
  268.      (cons
  269.       '(mlist)
  270.       (cons (pdis* (car p))
  271.         (mapcar 'rdis* (cdr p))))))
  272.  
  273. (defun psplit-field (p)                    ;modresult?
  274.   (let ((l (mapcar #'(lambda (q) (psplit-field1 (cdr q)))
  275.            (good-form (pfactor p))))    ;don't normalize lcfs?
  276.     ($algebraic t))
  277.        (if (null (cdr l)) (car l)
  278.        (do ((l l (cdr l))
  279.         (prim) (zeroes) (temp))
  280.            ((null l) (cons (or prim '|$splits in q|) zeroes))
  281.            (cond ((eq (caar l) 'linear)
  282.               (setq zeroes (cons (cdar l) zeroes)))
  283.              ((null prim)
  284.               (setq prim (caar l)
  285.                 zeroes (nconc (cdar l) zeroes)))
  286.              ((setq temp
  287.                 (primelmt (cons (car p) (cdr prim))
  288.                       (cons (car p) (cdaar l))
  289.                       *algvar)
  290.                 zeroes
  291.                 (nconc
  292.                  (mapcar
  293.                   #'(lambda (q)
  294.                     (ratgsubst (cadddr temp) (caaar l) q))
  295.                   (cdar l))
  296.                  (mapcar
  297.                   #'(lambda (q)
  298.                     (ratgsubst (caddr temp) (car prim) q))
  299.                   zeroes))
  300.                 prim (car temp))))))))
  301.  
  302. (defun plsolve (p)
  303.        (ratreduce (pterm (cdr p) 0)
  304.           (pminus (pterm (cdr p) 1))))
  305.  
  306.  
  307. (defun psplit-field1 (p)
  308.        ;;returns minimal poly and list of roots
  309.        ;;p must be square free
  310.    (*bind* ((minpoly (cons *algvar (cdr p)))
  311.        (zeroes) ($algebraic t)
  312.        ($ratalgdenom t))
  313.      (if (equal (cadr p) 1) (return (cons 'linear (plsolve p))))
  314.      (do ((polys (list p) )
  315.       (nminpoly)
  316.       (*nosplitf nil nil)
  317.       (alpha (cons (make-poly (car minpoly)) 1)))
  318.      ((null polys)
  319.       (cons minpoly zeroes))
  320.      (push alpha zeroes)
  321.      (putprop (car minpoly) (cdr minpoly) 'tellrat)
  322.      (rplaca polys
  323.          (car
  324.           (rquotient (pctimes (cdr alpha) (car polys))
  325.                  (pdifference
  326.                   (pctimes (cdr alpha) (pget (caar polys)))
  327.                   (car alpha)))))
  328.      (setq polys
  329.            (mapcan
  330.         #'(lambda (q) 
  331.               (cond ((equal (cadr q) 1)    ;;linear factor
  332.                  (push (plsolve q) zeroes)
  333.                  nil)            ;;flush linear factor
  334.                 ((list q))))
  335.         (mapcan #'(lambda (q)
  336.                   (let ((*ans) (*num 1) (*denom 1))
  337.                        (nreverse (pfactoralg1 q minpoly 0))))
  338.             polys)))
  339.      (when *nosplitf
  340.            (setq nminpoly (car *nosplitf)
  341.              *nosplitf (cdr *nosplitf))
  342.            (putprop *algvar (cdr nminpoly) 'tellrat)
  343.            (let ((beta
  344.               (plsolve (pgcd (cons (caar *nosplitf) (cdr minpoly))
  345.                      (exchangevar (car *nosplitf) *algvar)))))
  346.             (setq alpha (ratplus (cons (make-poly *algvar) 1)
  347.                      (rattimes (cons (f- (cadr *nosplitf)) 1)
  348.                            beta t)))
  349.             (setq zeroes
  350.               (mapcar
  351.                #'(lambda (q) (ratgsubst beta (car minpoly) q))
  352.                zeroes))
  353.             (setq polys
  354.               (mapcar
  355.                #'(lambda (q) (car (rgsubst beta (car minpoly) q)))
  356.                polys))
  357.             (setq minpoly
  358.               (cons *algvar (cdr nminpoly))))))))
  359.  
  360. (defun exchangevar (poly var)
  361.        (let ((newvar (gensym))
  362.          (ovar (car poly)))
  363.      (set newvar (f1+ (eval ovar)))
  364.      (pvsubst ovar newvar
  365.           (pvsubst var ovar
  366.                (pvsubst newvar var poly)))))
  367.  
  368. (defun rgsubst (val var p)    ;;generalized psubst substitutes any
  369.        (cond ((pcoefp p)
  370.           (cons p 1))    ;;expression for any var in p
  371.          ((eq var (car p))
  372.           (cond ((pzerop val)
  373.              (cons (pterm (cdr p) 0) 1))
  374.             ((do ((ld (cadr p) (car a))
  375.               (a (cdddr p) (cddr a))
  376.               (ans (cons (caddr p) 1)
  377.                    (ratplus
  378.                 (rattimes ans
  379.                       (ratexpt val
  380.                            (f- ld (car a)))
  381.                       t)
  382.                 (cons (cadr a) 1))))
  383.              ((null a) (rattimes ans (ratexpt val ld) t))))))
  384.          ((pointergp var (car p)) (cons p 1))
  385.          (t (let ((newsym (gensym)))
  386.           (set newsym (f1+ (symbol-value (car p))))
  387.           (rgsubst val newsym (pvsubst newsym var p))))))
  388.  
  389. (defun ratgsubst (val var rat) 
  390.        (ratquotient (rgsubst val var (car rat))
  391.             (rgsubst val var (cdr rat))))
  392.  
  393. (defun algnorm (f p)
  394.        (presult f p (car p)))
  395.  
  396. (DEFMFUN $algnorm (r p var)
  397.      (let ((varlist (list var))
  398.            (genvar))
  399.           (setq r (ratf r)
  400.             p (cadr (ratf p)))
  401.           (rdis* (cons (algnorm (cadr r) p)
  402.                (algnorm (cddr r) p)))))
  403.  
  404. (defun sqfrnorm (f p fvar)    ;;f must be sqfr, p is minpoly, fvar # pvar
  405.        (*bind* ((pvar (car p)))
  406.          (setq f (cdr (ordervar pvar (list f p)))  ;;new main var will be car of p
  407.            p (cadr f) f (car f))     ;make mainvar of f = mainvar(p)
  408.          (do ((i 0 (f1+ i))
  409.           (dif (pdifference (make-poly fvar) (make-poly (car p))))
  410.           (f f (pgsubst dif fvar f))
  411.           (res))
  412.          ((and (eq (car f) (car p))
  413.                (setq res (primpart (algnorm f p)))
  414.                (psqfrp res fvar))
  415.           (list res
  416.             (*bind* (($algebraic t))        ;;;modified f
  417.                   (putprop pvar (cdr p) 'tellrat)
  418.                   (pvsubst pvar (car p) f))
  419.             (car p)
  420.             p
  421.             i)))))
  422.  
  423. (defun primelmt (a b gvar &aux ($algebraic nil))
  424.        ;;a is a poly with coeff's in k(b)
  425.        ;;gvar is new variable
  426.        (let ((norm (sqfrnorm (cons gvar (cdr a)) b gvar))
  427.          (alpha) (beta) ($ratalgdenom t))
  428.         (rplaca norm (primpart (car norm)))
  429.         (putprop gvar (cdar norm) 'tellrat)
  430.         (setq $algebraic t
  431.           beta (pgcd (cadddr norm)
  432.                  (pvsubst (caddr norm)
  433.                       (car b)
  434.                       (cadr norm))))
  435.         (setq beta (plsolve beta)
  436.           alpha (ratplus (cons (make-poly gvar) 1)
  437.                  (rattimes (cons (f- (cadddr (cdr norm))) 1)
  438.                        beta t)))
  439.         (list (car norm)            ;;minimal poly
  440.           (pplus (make-poly a)    ;;new prim elm in old guys
  441.              (list (car b) 1 (f- (cadddr (cdr norm)))))
  442.           alpha beta)))            ;;in terms of gamma
  443.  
  444. (comment Discriminant of a Basis)
  445. (DEFMFUN $bdiscr n
  446.        (let ((varlist) (genvar))
  447.      (xcons (bdiscr (mapcar #'rform (listify (f1- n)))
  448.             (car (rform (arg n))))
  449.         (list 'mrat 'simp varlist genvar))))
  450.  
  451. (defun bdiscr (l minp)
  452.        (det
  453.     (mapcar #'(lambda (q1)
  454.             (mapcar #'(lambda (q2) (algtrace (ptimes q1 q2) minp)) l))
  455.         l)))
  456.  
  457. ;(declare-top(unspecial genvar varlist *var *num *ans alpha adn*))
  458.